home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / binomial.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  113 lines
  1. ;$Id: binomial.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       BINOMIAL
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the probabilty (bp) such that:
  11. ;                   Probability(X => v) = bp 
  12. ;       where X is a random variable from the cumulative binomial distribution
  13. ;       (Bernouli distribution).
  14. ;
  15. ; CATEGORY:
  16. ;       Statistics.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = Binomial(V, N, P)
  20. ;
  21. ; INPUTS:
  22. ;       V:    A non-negative integer specifying the minimal number of 
  23. ;             times an event E occurs in (N) independent performances.
  24. ;
  25. ;       N:    A non-negative integer specifying the number of performances.
  26. ;             If the number of performances exceeds 25, the Gaussian 
  27. ;             distribution is used to approximate the cumulative binomial 
  28. ;             distribution.
  29. ;
  30. ;       P:    A non-negative scalar, in the interval [0.0, 1.0],  of type 
  31. ;             float or double that specifies the probability of occurance 
  32. ;             or success of a single independent performance.
  33. ;
  34. ; EXAMPLES:
  35. ;       Compute the probability of obtaining at least two 6s in rolling a
  36. ;       die four times. The result should be 0.131944
  37. ;         result = binomial(2, 4, 1./6.)
  38. ;
  39. ;       Compute the probability of obtaining exactly two 6s in rolling a
  40. ;       die four times. The result should be 0.115741
  41. ;         result = binomial(2, 4, 1./6.) - binomial(3, 4, 1./6.)
  42. ;
  43. ;       Compute the probability of obtaining three or fewer 6s in rolling
  44. ;       a die four times. The result should be 0.999228
  45. ;         result = (binomial(0, 4, 1./6.) - binomial(1, 4, 1./6.)) + $
  46. ;                  (binomial(1, 4, 1./6.) - binomial(2, 4, 1./6.)) + $
  47. ;                  (binomial(2, 4, 1./6.) - binomial(3, 4, 1./6.)) + $
  48. ;                  (binomial(3, 4, 1./6.) - binomial(4, 4, 1./6.))
  49. ;
  50. ; PROCEDURE:
  51. ;       BINOMIAL computes the probability that an event E occurs at least
  52. ;       (V) times in (N) independent performances. The event E is assumed
  53. ;       to have a probability of occurance or success (P) in a single 
  54. ;       performance.
  55. ;
  56. ; REFERENCE:
  57. ;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
  58. ;       Erwin Kreyszig
  59. ;       ISBN 0-471-55380-8
  60. ;
  61. ; MODIFICATION HISTORY:
  62. ;       Modified by:  GGS, RSI, July 1994
  63. ;                     Minor changes to code. Rewrote documentation header.
  64. ;-
  65.  
  66. function N_BANG, n, min, fac1
  67.   ;If min and fac1 are undefined, then N_BANG returns n!.
  68.   ;Otherwise, fac1 * min * (min+1)....n is returned.
  69.   if n_elements(min) eq 0 then min = 2
  70.   if n_elements(fac1) eq 0 then fac = 1. $
  71.     else fac = fac1
  72.   if min gt n then return, fac
  73.   n1 =  n < 10
  74.   if min lt 11 then  $
  75.     for i = min, n1 do fac = i*fac
  76.   if (n lt 11.) then return, fac
  77.   n1 = 11 > min
  78.   ;Use logarithms to preserve precision.
  79.   return, fac * exp(total(alog(findgen(n-n1+1) + n1)))
  80. end
  81.  
  82. function binomial, v, n, p
  83.  
  84.   on_error, 2  ;Return to caller if error occurs.
  85.  
  86.   if p lt 0. or p gt 1. then message, $
  87.     'p must be in the interval [0.0, 1.0]'
  88.  
  89.   if v lt 0 then message, $
  90.     'v must be nonnegative.'
  91.  
  92.   if n lt 0 then message, $
  93.     'n must be nonnegative.'
  94.  
  95.   if v eq 0 then return,  1.0 $
  96.   else if n gt 25 then return, $
  97.                   1.0 - gauss_pdf((v-n*p)/sqrt(n*p*(1-p))) $
  98.   else if v gt n then return,  0.0
  99.  
  100.   nn = fix(n)
  101.   vv = fix(v)
  102.   n2 = vv < (nn - vv)
  103.   n3 = vv > (nn - vv)
  104.   n1 = N_BANG(nn, n3+1, p)
  105.   n1 =  n1/N_BANG(n2)
  106.   sum = n1 * p^(vv-1) * (1-p)^(nn-vv)
  107.   for i = vv+1, nn do begin
  108.     n1 = (nn-i+1) * n1/float(i)
  109.     sum = sum + n1 * p^(i-1) * (1-p)^(nn-i)
  110.   endfor
  111.   return, sum
  112. end
  113.